Categorical Attributes

• workclass - Which business sector do they work in?
• education - What level of education received?
• marital_status - What is their marriage history
• occupation - What do they do for a living
• relationship - Family member relation
• race - What is the subjects race
• gender - What is the subjects gender
• native_country - Where is the subject originally from
• income_bracket - Do they make over or under 50k/year

Continuous Attributes

• age - How old is the subject?
• fnlwgt - Sampling weight of observation
• education_num - numerical encoding of education variable
• capital_gain - income from investment sources, separate from wages/salary
• capital_loss - losses from investment sources, separate from wages/salary
• hours_per_week - How many hours a week did they work?

Data Quality Inspection

The first thing we must do is check and see if there are any NAs in our dataset, just to make sure to not mess up our analysis.

NA_sum <- sort(sapply(data, function(x) sum(is.na(x))), decreasing = TRUE)
data.frame((NA_sum))

Data Cleaning

Looks like we are doing ok here. The next issue we have in the dataset, is because of the way the csv was stored, some of the levels in our factors include leading and trailing whitespace. This is highly undesirable, so we must clean it up:

GetFactors <- function(df) {
    return(names(Filter(is.factor, df)))
}
FixLevels <- function(x) {
    levels(x) <- trimws(levels(x))
    return(x)
}
data[GetFactors(data)] <- lapply(data[GetFactors(data)], FixLevels)
pander(lapply(data[GetFactors(data)], levels))

• workclass: ?, Federal-gov, Local-gov, Never-worked, Private, Self-emp-inc, Self-emp-not-inc, State-gov and Without-pay
• education: 10th, 11th, 12th, 1st-4th, 5th-6th, 7th-8th, 9th, Assoc-acdm, Assoc-voc, Bachelors, Doctorate, HS-grad, Masters, Preschool, Prof-school and Some-college
• marital_status: Divorced, Married-AF-spouse, Married-civ-spouse, Married-spouse-absent, Never-married, Separated and Widowed
• occupation: ?, Adm-clerical, Armed-Forces, Craft-repair, Exec-managerial, Farming-fishing, Handlers-cleaners, Machine-op-inspct, Other-service, Priv-house-serv, Prof-specialty, Protective-serv, Sales, Tech-support and Transport-moving
• relationship: Husband, Not-in-family, Other-relative, Own-child, Unmarried and Wife
• race: Amer-Indian-Eskimo, Asian-Pac-Islander, Black, Other and White
• gender: Female and Male
• native_country: ?, Cambodia, Canada, China, Columbia, Cuba, Dominican-Republic, Ecuador, El-Salvador, England, France, Germany, Greece, Guatemala, Haiti, Holand-Netherlands, Honduras, Hong, Hungary, India, Iran, Ireland, Italy, Jamaica, Japan, Laos, Mexico, Nicaragua, Outlying-US(Guam-USVI-etc), Peru, Philippines, Poland, Portugal, Puerto-Rico, Scotland, South, Taiwan, Thailand, Trinadad&Tobago, United-States, Vietnam and Yugoslavia
• income_bracket: <=50K and >50K

Next, we need to reencode our data as factors. First, lets encode the education levels into factors with larger groups (for example 1st-12th grade should be no diploma, not a bunch of levels).

data$education <- fct_collapse(data$education, `No Diploma` = c("1st-4th", "5th-6th", 
    "7th-8th", "9th", "10th", "11th", "12th", "Preschool"), Associates = c("Assoc-acdm", 
    "Assoc-voc"), Diploma = c("Some-college", "HS-grad"))

Then the the income brackets:

data$income_bracket <- fct_collapse(data$income_bracket, small = "<=50K", large = ">50K")

Next, lets change the ? levels to something more useful:

levels(data$workclass)[levels(data$workclass) == "?"] <- "Other"
levels(data$occupation)[levels(data$occupation) == "?"] <- "Other-service"
levels(data$native_country)[levels(data$native_country) == "?"] <- "Other"
data[GetFactors(data)] <- lapply(data[GetFactors(data)], FixLevels)
pander(lapply(data[GetFactors(data)], levels))

• workclass: Other, Federal-gov, Local-gov, Never-worked, Private, Self-emp-inc, Self-emp-not-inc, State-gov and Without-pay
• education: No Diploma, Associates, Bachelors, Doctorate, Diploma, Masters and Prof-school
• marital_status: Divorced, Married-AF-spouse, Married-civ-spouse, Married-spouse-absent, Never-married, Separated and Widowed
• occupation: Other-service, Adm-clerical, Armed-Forces, Craft-repair, Exec-managerial, Farming-fishing, Handlers-cleaners, Machine-op-inspct, Priv-house-serv, Prof-specialty, Protective-serv, Sales, Tech-support and Transport-moving
• relationship: Husband, Not-in-family, Other-relative, Own-child, Unmarried and Wife
• race: Amer-Indian-Eskimo, Asian-Pac-Islander, Black, Other and White
• gender: Female and Male
• native_country: Other, Cambodia, Canada, China, Columbia, Cuba, Dominican-Republic, Ecuador, El-Salvador, England, France, Germany, Greece, Guatemala, Haiti, Holand-Netherlands, Honduras, Hong, Hungary, India, Iran, Ireland, Italy, Jamaica, Japan, Laos, Mexico, Nicaragua, Outlying-US(Guam-USVI-etc), Peru, Philippines, Poland, Portugal, Puerto-Rico, Scotland, South, Taiwan, Thailand, Trinadad&Tobago, United-States, Vietnam and Yugoslavia
• income_bracket: small and large

Next, lets remove the fnlwgt, education number, and capital gain and loss columns, as they are unneeded. We also need to rename some columns to be easier for us, and use the cut function to factorize our numeric variables

data <- data[, -c(3, 5, 11:12)]
data$age <- cut(data$age, breaks = c(15, 25, 45, 65, 100), labels = c("Young", 
    "Middleaged", "Senior", "Retired"))
data$hours_per_week <- cut(data$hours_per_week, breaks = c(0, 20, 40, 60, 80), 
    labels = c("part-time", "full-time", "hard-working", "need-a-life"))
library(magrittr)
data %>% skim_to_list() %$% factor %>% select(-c("ordered", "missing", "complete")) %>% 
    kable

Lets see the results:

Rule quality and inspection

Next, let us inspect the rules, and examine their quality

(quality(rulesorted))

inspectDT(rulesorted, caption = "Association Rules")

Plots

plot(rulesorted, method = "scatterplot", measure = c("support", "lift"), shading = "confidence", 
    engine = "htmlwidget")

plot(rulesorted, method = "graph", measure = "confidence", shading = "lift", 
    engine = "htmlwidget")

Parallel plot

plot(rulesorted, method = "paracoord", measure = "confidence", shading = "lift", 
    control = list(reorder = T))

Parallel Plot

Discussion

The width of the arrows represents support and the color represents confidence. We can see a few interesting things here. First of all, we have a striking level of confidence that being married gets you a large income, and then also interestingly we see that working as a private worker/housekeeper/server is a part of a lot of the rules related to being in the low income bracket (You can see this by looking where all the rules converge on LHS order 1).

plot(rulesorted, method = "two-key plot", measure = "lift", shading = "confidence", 
    engine = "htmlwidget")

Inspection

head(quality(rulesorted2))

inspectDT(rulesorted2, caption = "Alternate Association Rules")

Plots

plot(rulesorted2, method = "scatterplot", measure = c("support", "lift"), shading = "confidence", 
    engine = "htmlwidget")

plot(rulesorted2, method = "graph", measure = "confidence", shading = "lift", 
    engine = "htmlwidget")

plot(rulesorted2, method = "two-key plot", measure = "confidence", shading = "lift", 
    engine = "htmlwidget")

Small Income

Now that we’ve scanned a couple of different variations of apriori and generated rules accordingly. Lets split up the targeting in the rhs section to only target one of the income brackets at a time. First we’ll start with low income.

rule_small <- apriori(data, parameter = list(minlen = 2, supp = 0.12, conf = 0.95), 
    appearance = list(rhs = c("income_bracket=small"), default = "lhs"), control = list(verbose = F))
length(rule_small)

#> [1] 45

redundant <- is.redundant(rule_small)
rulep <- rule_small[redundant == FALSE]
rulesorted_small <- sort(rulep, by = "lift", decreasing = TRUE)
length(rulesorted_small)

#> [1] 25

head(quality(rulesorted_small))

inspectDT(rulesorted_small)

It looks like we’ve generated 25 rules for the low income bracket. Lets take a look at the scatter plot and find our releveant itemsets/rules

plot(rulesorted_small, method = "scatterplot", measure = c("confidence", "support"), 
    shading = "lift", engine = "htmlwidget")

We will discuss the results from small income targeting in the conclusion section.

Large Income

Our next step is to now target the large income bracket and see what rules shake out

rule_large <- apriori(data, parameter = list(minlen = 2, supp = 0.05, conf = 0.55), 
    appearance = list(rhs = c("income_bracket=large"), default = "lhs"), control = list(verbose = F))
length(rule_large)

#> [1] 53

redundant <- is.redundant(rule_large)
rulep <- rule_large[redundant == FALSE]
rulesorted_large <- sort(rulep, by = "lift", decreasing = TRUE)
length(rulesorted_large)

#> [1] 15

head(quality(rulesorted_large))

inspectDT(rulesorted_large)

It looks like we’ve generated 25 rules for the low income bracket. Lets take a look at the scatter plot and find our releveant itemsets/rules

plot(rulesorted_large, method = "scatterplot", measure = c("confidence", "support"), 
    shading = "lift", engine = "htmlwidget")

We will discuss the results from large income targeting in the conclusion section.

Small Income

Interesting small Income rules:

• Not being married is associated with a low income. Out of our 25 best rules, 16 of them involved not being married, and our five highest lift rules all involved never married. This is clearly not a coincidence, however what this means and the cause of this is not so clear. We propose a few hypotheses about this:

  • Married people, motivated by love, work harder and make more money

  • Shallow society deems people without money harder to marry

• People who are married got married because they are more social, and thus they are given more raises, perform better in interviews, etc. In contrast, people who were not sociable enough to meet someone may not perform as well in interviews etc., giving them fewer opportunities to get a job

• Having just a diploma from high school does not get you out of the lower income bracket. This is unsurprising, as today you need in general a minimum of a college education to start making good money (and even that is difficult). 5 out of our 25 rules involved highest education level = high school diploma

• 12 out of 25 rules involved being young. This makes perfect sense, as entry level jobs which young people get do not pay very well.

• 8 out of 25 rules involved being in the private /housekeeping/server field. This makes sense, as these workers are severly underpaid

• 4 out of 25 rules involved having a child. If you have a child, you should be spending more time with them, and maybe focusing a little less on your career. Therefore, you will probably earn less. This pattern makes sense.

• The rules for our small income target have a support of 12% and confidence of 95%. This means that (never-married, Diploma, Young, Private/housekeeper/server, and child) are observed in 12% of the dataset. We therefore believe that 95% of the time any of the attributes (never-married, Diploma, Young, Private/housekeeper/server, and child) are present, it will result in a small income.

Large Income

Interesting large Income rules:

• Thankfully, opposite to the small income analysis, large income targeting shows that being married greatly influences your income as it shows up in 9 of the 15 rules present. Meaning again that by working as a team, you get more out of life. So go find yourself a husband/wife.

• The next highest frequent itemset is hard-working that shows up in 8 of the 15 rules. Which makes perfect sense that the more time you put into a job, the likelier you are for promotion and advancement within an organization.

• Education = Bachelors , native_country = United-States, relationship = husband are present in 5 out of the 15 rules.

• The rules for our large income target have a support of 05% and confidence of 55%. This means that (married, hard-working, Bachelors, United-States, Husband) are observed in 5% of the dataset. We therefore believe that 55% of the time any of the attributes (married, hard-working, Bachelors, United-States, Husband) are present, it will result in a large income.

Building a Grid

There are two parameters to optimize in the apriori algorithm: confidence and support. Normally, these are chosen using an in depth knowledge of the data combined with an understanding of the goals of the study. Because this research is purely academic, and we are new to the apriori algorithm, we do not have a prior understanding to tell us what support and confidence will be the most useful for us. So, we must optimize both. We will begin by creating two objects, conf and supp, for confidence and support respectively. Each will be vectors from 0.01 to 0.99:

x <- 1:99
conf <- supp <- x/100

Next, we will create a data frame containing all possible combinations of support and confidence:

init <- data.frame(conf = conf, supp = supp)
params <- expand.grid(init)
nrow(params)

#> [1] 9801

Defining the Search Function

Next, we would like to create a function which we can iterate through. We would like the function to take in a value of support and confidence, maybe a value of appearance (but with defaults for our specific case) and run it through the apriori algorithm. Then, we would like it to prune redundant results from the rule, get the length and average lift of the ruleset, and return those in a convenient labeled vector:

search <- function(sup, con, appearance = list(rhs = c("income_bracket=small", 
    "\n                                            income_bracket=large"), default = "lhs")) {
    rules <- apriori(transactions, parameter = list(minlen = 2, supp = sup, 
        conf = con), control = list(verbose = f))
    redundant <- is.redundant(rules)
    rulepruned <- rules[redundant == false]
    rulelength <- length(rules)
    return(c(support = sup, conf = con, nrules = rulelength, avlift = mean(quality(rulepruned)$lift)))
}

Performing the grid search

We will perform the grid search in parallel, in two ways: using R’s built in parallel library, utilizing a FORK cluster (not available on Windows), and using doParallel and foreach, with a socket cluster (cross platform).

library(parallel)
gs <- mcmapply(search, params[[1]], params[[2]], mc.cores = 12L)
results <- data.frame(t(gs))

libary(doParallel)
library(foreach)
numCores <- parallel::detectCores()
workers <- makeCluster(numCores - 1)
invisible(clusterEvalQ(workers, library(arules)))
registerDoParallel(workers)

results <- foreach(i = 1:nrow(params), .combine = rbind) %dopar% {
    search(params[i, 1], params[i, 2])
}

Viewing the Results

We now must define what it means to be digestible. In my opinion, between 20 and 25 rules is clearly digestible to anyone. So we will find the rules created that give us that, and we want to also find the rulesets with the highest lift. We will do so using dplyr:

digest <- results %>% filter(nrules <= 25 & nrules >= 20) %>% filter(conf > 
    0.35) %>% filter(support > 0.35) %>% arrange(desc(avlift, nrules, support, 
    conf))
digest

Looks like we have a clear winner. We will now try